{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from [this Jupyter notebook](https://nbviewer.jupyter.org/github/krischer/seismo_live/blob/master/notebooks/Computational%20Seismology/The Finite-Difference Method/fd_taylor_operators_advanced.ipynb) by Heiner Igel ([@heinerigel](https://github.com/heinerigel)), Lion Krischer ([@krischer](https://github.com/krischer)) which is a supplemenatry material to the book [Computational Seismology: A Practical Introduction](http://www.computational-seismology.org/),  additional modifications by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Improving the accuracy of FD operators based on Taylor series expansion\n",
    "\n",
    "After estimating the accuracy of FD operators by calculating the remainders of Taylor series expansions, the next obvious question is how could we improve the accuracy of FD approximations.\n",
    "\n",
    "## Estimation of FD operators by Taylor series expansion\n",
    "\n",
    "In [this lession](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/03_Intro_finite_differences/1_fd_intro.ipynb), we introduced the 3-point FD approximation to the second derivative of a function:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx \\frac{f(x+dx)-2f(x)+f(x-dx)}{dx^2} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "from a Taylor series expansion. Here, I want to introduce a more systematic approach to derive this operator, which is required for a later improvement of the FD operator accuracy.\n",
    "\n",
    "The Taylor expansion of $f(x + dx)$ around $x$ is defined as \n",
    "\n",
    "$$\n",
    "f(x+dx)=\\sum_{n=0}^\\infty \\frac{f^{(n)}(x)}{n!}dx^{n}  \n",
    "$$\n",
    "\n",
    "Finite-difference operators can be calculated by seeking weights (here: $a$, $b$, $c$) with which function values have to be multiplied to obtain an interpolation or a derivative. \n",
    "\n",
    "### Example: second derivative\n",
    "\n",
    "We want to approximate the second derivative of a function based on function values at the points $f(x+dx)$, $f(x)$ and $f(x-dx)$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx a ~ f(x+dx) + b ~ f(x) + c ~ f(x-dx)\n",
    "\\end{equation}\n",
    "\n",
    "The individual terms in the derivative can be expressed by the following Taylor series expansion up to the quadratic term:\n",
    "\n",
    "\\begin{align}\n",
    "a ~ f(x + dx) & \\ = \\ a ~ \\left[ ~  f(x) +  \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{1}{2!}  \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2   + \\dotsc  ~ \\right] \\nonumber\\\\\n",
    "b ~ f(x) & \\  = \\ b ~ \\left[ ~  f(x)  ~ \\right] \\nonumber\\\\\n",
    "c ~ f(x - dx) & \\ = \\ c ~ \\left[ ~  f(x) -  \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{1}{2!}  \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2   - \\dotsc  ~ \\right] \\nonumber\n",
    "\\end{align}\n",
    "\n",
    "Inserting these approximations into eq. (1) and sorting of the terms leads to \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx a ~ f(x+dx) + b ~ f(x) + c ~ f(x-dx) \\approx \\nonumber\n",
    "\\end{equation}\n",
    "\\begin{align}\n",
    "f(x) ~ \\biggl[&a  ~~+~~            ~~~~b           &+~~  c\\biggr] \\notag\\\\\n",
    "+ ~ dx ~ \\frac{\\partial f(x)}{\\partial x} ~ \\biggl[&a  ~~\\phantom{+}~~ \\phantom{b}  &-~~  c\\biggr] \\notag\\\\\n",
    "+ ~ \\frac{dx^2}{2!} ~ \\frac{\\partial^2 f(x)}{\\partial x^2} ~ \\biggl[&a  ~~\\phantom{+}~~ \\phantom{b}  &+~~  c\\biggr] \\notag\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "The RHS of this equation can only be an approximation to $\\frac{\\partial^2 f(x)}{\\partial x^2}$, if the following conditions are fulfilled:\n",
    "\n",
    "\\begin{align}\n",
    "&a  ~~+~~            ~~~~b           &+~~  c & = & 0 \\notag\\\\\n",
    "&a  ~~\\phantom{+}~~ \\phantom{b}  &-~~  c & = & 0 \\notag\\\\\n",
    "&a  ~~\\phantom{+}~~ \\phantom{b}  &+~~  c & = & \\frac{2!}{\\mathrm{d}x^2}\\notag\n",
    "\\end{align}\n",
    "\n",
    "This system of linear equations can be expressed in matrix form as:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{pmatrix}\n",
    "  1 & 1 & 1 \\\\\n",
    "  1 & 0 & -1 \\\\\n",
    "  1 & 0 & 1\n",
    " \\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "  a\\\\\n",
    " b \\\\\n",
    " c\n",
    " \\end{pmatrix}\n",
    " =\n",
    " \\begin{pmatrix}\n",
    "  0\\\\\n",
    " 0 \\\\\n",
    " \\frac{2!}{dx^2}\n",
    " \\end{pmatrix}\n",
    " \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "or \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{Ax ~= ~s}\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "with the system matrix $\\mathbf{A}$, the vector $\\mathbf{x}$ containing the desired FD weights $a,\\; b,\\; c$ and the solution vector $\\mathbf{s}$. Formally, the FD operator weights can be estimated by matrix inversion:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{x ~= A^{-1}~s}\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Let's try to estimate the FD weighting coefficients by symbolic matrix inversion with `SymPy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import SymPy libraries\n",
    "from sympy import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define symbols\n",
    "a, b, c, dx = symbols('a b c dx')\n",
    "\n",
    "# In SymPy we define the matrix equation by a matrix consisting of the \n",
    "# system matrix A and an additional column defining the solution vector s\n",
    "\n",
    "\n",
    "# Symbolic solution of the linear system can be easily obtained by linsolve\n",
    "linsolve(system, (a, b, c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the resulting weighting coefficients are:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{pmatrix}\n",
    "a \\\\\n",
    "b\\\\\n",
    "c\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "\\frac{1}{\\mathrm{d}x^2} \\\\\n",
    " - \\frac{2}{\\mathrm{d}x^2} \\\\\n",
    "  \\frac{1}{\\mathrm{d}x^2}\n",
    "\\end{pmatrix}\n",
    "\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Inserting the coefficients $a,\\; b,\\; c$ in eq. (1) leads to \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx \\frac{f(x+dx)-2f(x)+f(x-dx)}{dx^2}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "which is the well known 3-point operator for the 2nd derivative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Improving the accuracy of the 2nd derivative FD approximation\n",
    "\n",
    "We can use the above approach to estimate higher order and therefore more accurate FD operators. Let's assume we want to approximate the second derivative by a 5- instead of a 3-point operator, at the points $f(x+2dx)$, $f(x+dx)$, $f(x)$, $f(x-dx)$ and $f(x-2dx)$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx a ~ f(x+2dx) + b ~ f(x+dx) + c ~ f(x) + d ~ f(x-dx) + e ~ f(x-2dx) \\notag\n",
    "\\end{equation}\n",
    "\n",
    "To solve for the 5 unknown FD weights, we have to expand 5 terms of the Taylor series, which leads to a higher error order and therefore a more accurate FD operator compared to the 3-point operator. The resulting Taylor series expansions are:\n",
    "\n",
    "\\begin{align}\n",
    "a ~ f(x + 2dx) & \\ = \\ a ~ \\left[ ~  f(x) +  2 \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{4}{2}  \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2 +  \\frac{8}{6}  \\frac{\\partial^3 f(x)}{\\partial x^3} dx^3 +  \\frac{16}{24}  \\frac{\\partial^4 f(x)}{\\partial x^4} dx^4   + \\dotsc  ~ \\right] \\nonumber\\\\\n",
    "b ~ f(x + dx) & \\ = \\ b ~ \\left[ ~  f(x) +  \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{1}{2} \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2 +  \\frac{1}{6} \\frac{\\partial^3 f(x)}{\\partial x^3} dx^3 +  \\frac{1}{24}  \\frac{\\partial^4 f(x)}{\\partial x^4} dx^4   + \\dotsc  ~ \\right] \\nonumber\\\\\n",
    "c ~ f(x) & \\  = \\ c ~ \\left[ ~  f(x)  ~ \\right] \\nonumber\\\\\n",
    "d ~ f(x - dx) & \\ = \\ d ~ \\left[ ~  f(x) -  \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{1}{2} \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2 -  \\frac{1}{6} \\frac{\\partial^3 f(x)}{\\partial x^3} dx^3 +  \\frac{1}{24}  \\frac{\\partial^4 f(x)}{\\partial x^4} dx^4   + \\dotsc  ~ \\right] \\nonumber\\\\\n",
    "e ~ f(x - 2dx) & \\ = \\ e ~ \\left[ ~  f(x) -  2 \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{4}{2}  \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2 -  \\frac{8}{6} \\frac{\\partial^3 f(x)}{\\partial x^3} dx^3 +  \\frac{16}{24} \\frac{\\partial^4 f(x)}{\\partial x^4} dx^4   + \\dotsc  ~ \\right] \\nonumber\n",
    "\\end{align}\n",
    "\n",
    "Similar to the 3-point operator, we insert the terms into the 5-point FD approximation and sort the terms\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx a ~ f(x+2dx) + b ~ f(x+dx) + c ~ f(x) + d ~ f(x-dx) + e ~ f(x-2dx)  \\approx \\nonumber\n",
    "\\end{equation}\n",
    "\\begin{align}\n",
    "f(x) ~ \\biggl[&a  ~~+~~   b   ~~+~~ c ~~+~~  d ~~+~~  e~~           \\biggr] \\notag\\\\\n",
    "+ ~ dx ~ \\frac{\\partial f(x)}{\\partial x} ~ \\biggl[&2a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~-~~  d ~~-~~  2e~~           \\biggr] \\notag\\\\\n",
    "+ ~ \\frac{dx^2}{2} ~ \\frac{\\partial^2 f(x)}{\\partial x^2} ~ \\biggl[&4a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~+~~  d ~~+~~  4e~~           \\biggr] \\notag\\\\\n",
    "+ ~ \\frac{dx^3}{6} ~ \\frac{\\partial^3 f(x)}{\\partial x^3} ~ \\biggl[&8a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~-~~  d ~~-~~  8e~~           \\biggr] \\notag\\\\\n",
    "+ ~ \\frac{dx^4}{24} ~ \\frac{\\partial^4 f(x)}{\\partial x^4} ~ \\biggl[&16a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~+~~  d ~~+~~  16e~~           \\biggr] \\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "leading to the following conditions:\n",
    "\n",
    "\\begin{align}\n",
    "\\biggl[&a  ~~+~~   b   ~~+~~ c ~~+~~  d ~~+~~  e~~           \\biggr] ~=~ 0 \\notag\\\\\n",
    "\\biggl[&2a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~-~~  d ~~-~~  2e~~           \\biggr] ~=~ 0 \\notag\\\\\n",
    "\\biggl[&4a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~+~~  d ~~+~~  4e~~           \\biggr] ~=~ \\frac{2}{dx^2} \\notag\\\\\n",
    "\\biggl[&8a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~-~~  d ~~-~~  8e~~           \\biggr] ~=~ 0\\notag\\\\\n",
    "\\biggl[&16a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~+~~  d ~~+~~  16e~~           \\biggr] ~=~ 0\\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "or as matrix equation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{pmatrix}\n",
    "  1 & 1 & 1 & 1 & 1 \\\\\n",
    "  2 & 1 & 0 & -1 & -2 \\\\\n",
    "  4 & 1 & 0 & 1 & 4 \\\\\n",
    "  8 & 1 & 0 & -1 & -8 \\\\\n",
    "  16 & 1 & 0 & 1 & 16 \\\\\n",
    " \\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "  a\\\\\n",
    " b \\\\\n",
    " c \\\\\n",
    " d \\\\\n",
    " e \\\\\n",
    " \\end{pmatrix}\n",
    " =\n",
    " \\begin{pmatrix}\n",
    "  0\\\\\n",
    " 0 \\\\\n",
    " \\frac{2!}{dx^2} \\\\\n",
    " 0 \\\\\n",
    " 0\n",
    " \\end{pmatrix}\n",
    " \\nonumber\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define symbols\n",
    "a, b, c, d, e, dx = symbols('a b c d e dx')\n",
    "\n",
    "# In SymPy we define the matrix equation by a matrix consisting of the \n",
    "# system matrix A and an additional column defining the solution vector s\n",
    "\n",
    "\n",
    "# Symbolic solution of the linear system can be easily obtained by linsolve\n",
    "linsolve(system, (a, b, c, d, e))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the resulting weighting coefficients for the 5-point FD operator are:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "a \\\\\n",
    "b\\\\\n",
    "c\\\\\n",
    "d\\\\\n",
    "e\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "-\\frac{1}{12\\mathrm{d}x^2} \\\\\n",
    " \\frac{4}{3\\mathrm{d}x^2} \\\\\n",
    "- \\frac{5}{2\\mathrm{d}x^2}\\\\\n",
    "\\frac{4}{3\\mathrm{d}x^2}\\\\\n",
    "-\\frac{1}{12\\mathrm{d}x^2}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "##### Exercise\n",
    "\n",
    "To investigate the effect of the 5-point operator on the finite-difference based solution of the 1D acoustic wave equation:\n",
    "\n",
    "* Implement the 5-point operator to approximate the spatial 2nd derivative in the 1D acoustic FD modelling code. Leave an option to also use the 3-point operator\n",
    "* Compare the impact of the 3-point and 5-point operators on the modelling result by plotting the analytical seismograms $d^{analyt}$ and FD seismograms $d^{FD}$ as well as their differences, respectively. Use $v_0 = 333\\; m/s$, $dx = 2\\; m$ and $dt = 1\\;ms$ as modelling parameters. Quantify the data misfit also in terms of the L2-norm of the data residuals:\n",
    "\\begin{equation}\n",
    "E = \\frac{1}{2}\\sum_{i=0}^{nt}\\biggl(d^{FD}_i-d^{analyt}_i\\biggr)^2\\nonumber\n",
    "\\end{equation}\n",
    "where $nt$ denotes the number of time samples\n",
    "* Benchmark and compare the runtimes of the 3-point and 5-point operator by using the magic function `%timeit` and the same spatio-temporal model discretizations $dx$ and $dt$ in order to achieve comparable results\n",
    "* Estimate the factor $\\zeta$ in the Courant criterion:\n",
    "\\begin{equation}\n",
    "dt \\le \\frac{dx}{\\zeta v_{max}}\\nonumber\n",
    "\\end{equation}\n",
    "for the 5-point operator by forward modelling. Start with the modelling parameters $v_{max} = 333\\; m/s$ and $dx = 0.5\\; m$ and the Courant criterion for the 3pt operator ($\\zeta=1$). Vary $\\zeta$ during different forward modelling runs, until you resolved the transition between stable and unstable FD solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Libraries \n",
    "# ----------------\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "xmax = 500 # maximum spatial extension of the 1D model (m)\n",
    "dx   = 0.5 # grid point distance in x-direction\n",
    "\n",
    "tmax = 1.001   # maximum recording time of the seismogram (s)\n",
    "dt   = 0.0010  # time step\n",
    "\n",
    "vp0  = 333.   # P-wave speed in medium (m/s)\n",
    "\n",
    "# acquisition geometry\n",
    "xr = 365.0 # receiver position (m)\n",
    "xsrc = 249.5 # source position (m)\n",
    "\n",
    "f0   = 25. # dominant frequency of the source (Hz)\n",
    "t0   = 4. / f0 # source time shift (s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1D Wave Propagation (Finite Difference Solution) \n",
    "# ------------------------------------------------\n",
    "def FD_1D_acoustic(dt,dx,op):\n",
    "        \n",
    "    nx = (int)(xmax/dx) # number of grid points in x-direction        \n",
    "    nt = (int)(tmax/dt) # maximum number of time steps            \n",
    "    \n",
    "    ir = (int)(xr/dx)      # receiver location in grid in x-direction    \n",
    "    isrc = (int)(xsrc/dx)  # source location in grid in x-direction\n",
    "\n",
    "    # Source time function (Gaussian)\n",
    "    # -------------------------------\n",
    "    src  = np.zeros(nt + 1)\n",
    "    time = np.linspace(0 * dt, nt * dt, nt)\n",
    "\n",
    "    # 1st derivative of a Gaussian\n",
    "    src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "\n",
    "    # Analytical solution\n",
    "    # -------------------\n",
    "    G    = time * 0.\n",
    "\n",
    "    # Initialize coordinates\n",
    "    # ----------------------\n",
    "    x    = np.arange(nx)\n",
    "    x    = x * dx       # coordinate in x-direction\n",
    "\n",
    "    for it in range(nt): # Calculate Green's function (Heaviside function)\n",
    "        if (time[it] - np.abs(x[ir] - x[isrc]) / vp0) >= 0:\n",
    "            G[it] = 1. / (2 * vp0)\n",
    "    seis_ana   = np.convolve(G, src * dt)\n",
    "    seis_ana   = seis_ana[0:nt]\n",
    "    \n",
    "    # Initialize empty pressure arrays\n",
    "    # --------------------------------\n",
    "    p    = np.zeros(nx) # p at time n (now)\n",
    "    pold = np.zeros(nx) # p at time n-1 (past)\n",
    "    pnew = np.zeros(nx) # p at time n+1 (present)\n",
    "    d2px = np.zeros(nx) # 2nd space derivative of p\n",
    "\n",
    "    # Initialize model (assume homogeneous model)\n",
    "    # -------------------------------------------\n",
    "    vp    = np.zeros(nx)\n",
    "    vp    = vp + vp0       # initialize wave velocity in model\n",
    "\n",
    "    # Initialize empty seismogram\n",
    "    # ---------------------------\n",
    "    seis = np.zeros(nt)    \n",
    "    \n",
    "    # Calculate Partial Derivatives\n",
    "    # -----------------------------\n",
    "    for it in range(nt):\n",
    "    \n",
    "        # FD approximation of spatial derivative by 3 point operator\n",
    "        if(op==3):\n",
    "            for i in range(1, nx - 1):\n",
    "                d2px[i] = (p[i + 1] - 2 * p[i] + p[i - 1]) / dx ** 2\n",
    "        \n",
    "        # ADD YOUR 5-POINT FD OPERATOR HERE \n",
    "        #if(op==5):\n",
    "                \n",
    "        # Time Extrapolation\n",
    "        # ------------------\n",
    "        pnew = 2 * p - pold + vp ** 2 * dt ** 2 * d2px\n",
    "\n",
    "        # Add Source Term at isrc\n",
    "        # -----------------------\n",
    "        # Absolute pressure w.r.t analytical solution\n",
    "        pnew[isrc] = pnew[isrc] + src[it] / dx * dt ** 2\n",
    "                \n",
    "        # Remap Time Levels\n",
    "        # -----------------\n",
    "        pold, p = p, pnew\n",
    "    \n",
    "        # Output of Seismogram\n",
    "        # --------------------\n",
    "        seis[it] = p[ir]\n",
    "        \n",
    "    return time, seis, seis_ana"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = 2.0 # grid point distance in x-direction (m)\n",
    "dt = 0.0010  # time step (s)\n",
    "\n",
    "# run FD code with 3-point operator\n",
    "op = 3 # use 3-point operator\n",
    "time, seis_3pt, seis_ana = FD_1D_acoustic(dt,dx,op)\n",
    "\n",
    "# run FD code with 5-point operator\n",
    "op = 5 # use 5-point operator\n",
    "# RUN FD CODE WITH 5-POINT OPERATOR HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare FD Seismogram with analytical solution\n",
    "# ---------------------------------------------- \n",
    "# Define figure size\n",
    "rcParams['figure.figsize'] = 12, 5\n",
    "\n",
    "# plot solution of 3-point FD operator\n",
    "plt.plot(time, seis_3pt, 'b-',lw=3,label=\"FD solution (3pt operator)\")\n",
    "\n",
    "# PLOT SOLUTION OF 5-POINT FD OPERATOR HERE!\n",
    "\n",
    "# plot analytical solution\n",
    "Analy_seis = plt.plot(time,seis_ana,'r--',lw=3,label=\"Analytical solution\") \n",
    "\n",
    "# plot difference between analytical and 3pt operator FD solution\n",
    "plt.plot(time, 10*(seis_3pt-seis_ana), 'g-',lw=3,label=\"3pt operator error (10x)\")\n",
    "\n",
    "# PLOT DIFFERENCE BETWEEN ANALYTICAL AND 5PT OPERATOR SOLUTION HERE!\n",
    "\n",
    "plt.xlim(time[0], time[-1])\n",
    "plt.title('Seismogram comparison')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# Calculate L2-norm of the data residuals between FD and analytical solution\n",
    "# 3-point operator\n",
    "E_3pt = np.sum((seis_3pt-seis_ana)**2)\n",
    "print(\"E_3pt = \", E_3pt)\n",
    "\n",
    "# ADD L2-NORM for 5-POINT OPERATOR HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = 2.0 # grid point distance in x-direction\n",
    "dt = 0.0010  # time step\n",
    "\n",
    "# Benchmark runtime of 3pt-FD code\n",
    "op = 3 # use 3-point operator\n",
    "%timeit time, seis_3pt, seis_ana = FD_1D_acoustic(dt,dx,op)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BENCHMARK RUNTIME OF 5PT-FD CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate Courant criterion for 5pt-FD operator\n",
    "# ----------------------------------------------\n",
    "\n",
    "dx = 0.5 # grid point distance in x-direction\n",
    "zeta = 1.0\n",
    "dt = dx / (zeta*vp0)\n",
    "op = 5 # use 5-point operator\n",
    "# RUN FD CODE WITH 5-POINT OPERATOR HERE!\n",
    "\n",
    "# PLOT SOLUTION OF 5-POINT FD OPERATOR HERE!\n",
    "\n",
    "# plot analytical solution\n",
    "Analy_seis = plt.plot(time,seis_ana,'r--',lw=3,label=\"Analytical solution\") \n",
    "\n",
    "plt.xlim(time[0], time[-1])\n",
    "plt.title('Seismogram comparison')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* Use the Taylor series to calculate finite-difference operators\n",
    "* How to derive high-order finite-difference operators\n",
    "* Derivation and application of spatial 5-point operator to approximate the spatial second derivative in the 1D acoustic wave equation"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
